*robustness checks

**ROBUSTNESS 1
** resp characteristics as controls
   global mod1 = "i.sex  i.religion  i.residence i.occupation i.education i.language i.friends i.team_support i.children  i.country##i.legal_status i.gender i.region_grouped_it i.age_new i.educ_new i.income_hml_it i.empl_stat_eu i.voting_18 [pweight=weight]"
   reg outcome $mod1 , baselevels $ses coeflegend	 
     estimate store reg
	 estimate restore reg
	 restrictamce
*same graph*
*graph

coefplot (reg, keep( *.sex  *.religion  *.residence *.occupation *.education *.language *.friends  *.team_support *.children 1.country) pstyle(p1)) ///  
(Romania, rename ( (1) = R)  \ ///
country2, rename ( (1) = C)  \ ///
country4, rename ( (1) = P) \ ///
country5, rename ( (1) = S) pstyle(p1)) ///
(reg, keep (1.legal_status ) pstyle(p1)) ///
(legal, rename ( (1) = L) pstyle(p1)), ///
headings(1.sex = "{bf: Child gender}" 1.residence = "{bf: Parents' length of residence}" 1.occupation = "{bf: Parents' occupation}"  ///
1.education = "{bf: Parents' education}" 1.religion = "{bf: Parents' religion}" 1.language = "{bf: Parents' Italian proficiency}" 1.friends = "{bf: Family friends Origin}" 1.children = "{bf:Total number of children}" 1.team_support = "{bf:Parents' team support}" 1.country = "{bf: Parents' country of origin}" 1.legal_status = "{bf:Parents' legal status}") ///
xline (-0.4 -0.3 -0.2 -0.1 0.1 0.2, lstyle(grid)) xline (0, lpattern(dash)) ciopts(recast(rcap)) baselevels ///
coeflabels(    R = "Romania"  ///
               C = "China"  ///
               P = "Pakistan" ///
			   S = "Senegal"  ///
			   L = "Regular", labsize(vsmall)) msymbol(s) mfcolor(white)  offset(0.05) /// 
			  graphregion(margin(zero)) xlabel( -0.2 "-20" -0.1 "-10" 0 0.1 "10" 0.2 "20", valuelabel) xlabel(,labsize(small)) /// 
			    legend(off)   xtitle("Change in probability (perc. points)") graphregion(color(white)) xtitle(, size(small)) 	///
				order(*.occupation 1.legal_status L   *.residence  *.religion  *.language *.education *.friends  *.team_support *.children 1.country R C P S *.sex)
				graph export "FigS3_robustness.png", replace
				graph export "FigS3_robustness.eps", replace
**-----------------------------------------------------------------------------------------------------------------------------------------------

*fixed effects
    xtreg outcome $mod ,  $ses i(caseid) fe  
    estimate store reg
    restrictamce
	
	coefplot (reg, keep( *.sex  *.religion  *.residence *.occupation *.education *.language *.friends  *.team_support *.children 1.country) pstyle(p1)) ///  
(Romania, rename ( (1) = R)  \ ///
country2, rename ( (1) = C)  \ ///
country4, rename ( (1) = P) \ ///
country5, rename ( (1) = S) pstyle(p1)) ///
(reg, keep (1.legal_status ) pstyle(p1)) ///
(legal, rename ( (1) = L) pstyle(p1)), ///
headings(1.sex = "{bf: Child gender}" 1.residence = "{bf: Parents' length of residence}" 1.occupation = "{bf: Parents' occupation}"  ///
1.education = "{bf: Parents' education}" 1.religion = "{bf: Parents' religion}" 1.language = "{bf: Parents' Italian proficiency}" 1.friends = "{bf: Family friends Origin}" 1.children = "{bf:Total number of children}" 1.team_support = "{bf:Parents' team support}" 1.country = "{bf: Parents' country of origin}" 1.legal_status = "{bf:Parents' legal status}") ///
xline (-0.4 -0.3 -0.2 -0.1 0.1 0.2, lstyle(grid)) xline (0, lpattern(dash)) ciopts(recast(rcap)) baselevels ///
coeflabels(    R = "Romania"  ///
               C = "China"  ///
               P = "Pakistan" ///
			   S = "Senegal"  ///
			   L = "Regular", labsize(vsmall)) msymbol(s) mfcolor(white)  offset(0.05) /// 
			  graphregion(margin(zero)) xlabel( -0.2 "-20" -0.1 "-10" 0 0.1 "10" 0.2 "20", valuelabel) xlabel(,labsize(small)) /// 
			    legend(off)   xtitle("Change in probability (perc. points)") graphregion(color(white)) xtitle(, size(small)) 	///
				order(*.occupation 1.legal_status L   *.residence  *.religion  *.language *.education *.friends  *.team_support *.children 1.country R C P S *.sex) subtitle(B)
				graph save fixed_effects, replace
				
**-----------------------------------------------------------------------------------------------------------------------------------------------
				
*random effects
      global modre = "i.sex  i.religion  i.residence i.occupation i.education i.language i.friends i.team_support i.children  i.country##i.legal_status"
      xtreg outcome $modre ,   i(caseid) re cl(caseid)
	  estimate store reg
      restrictamce
	  
	  coefplot (reg, keep( *.sex  *.religion  *.residence *.occupation *.education *.language *.friends  *.team_support *.children 1.country) pstyle(p1)) ///  
(Romania, rename ( (1) = R)  \ ///
country2, rename ( (1) = C)  \ ///
country4, rename ( (1) = P) \ ///
country5, rename ( (1) = S) pstyle(p1)) ///
(reg, keep (1.legal_status ) pstyle(p1)) ///
(legal, rename ( (1) = L) pstyle(p1)), ///
headings(1.sex = "{bf: Child gender}" 1.residence = "{bf: Parents' length of residence}" 1.occupation = "{bf: Parents' occupation}"  ///
1.education = "{bf: Parents' education}" 1.religion = "{bf: Parents' religion}" 1.language = "{bf: Parents' Italian proficiency}" 1.friends = "{bf: Family friends Origin}" 1.children = "{bf:Total number of children}" 1.team_support = "{bf:Parents' team support}" 1.country = "{bf: Parents' country of origin}" 1.legal_status = "{bf:Parents' legal status}") ///
xline (-0.4 -0.3 -0.2 -0.1 0.1 0.2, lstyle(grid)) xline (0, lpattern(dash)) ciopts(recast(rcap)) baselevels ///
coeflabels(    R = "Romania"  ///
               C = "China"  ///
               P = "Pakistan" ///
			   S = "Senegal"  ///
			   L = "Regular", labsize(vsmall)) msymbol(s) mfcolor(white)  offset(0.05) /// 
			  graphregion(margin(zero)) xlabel( -0.2 "-20" -0.1 "-10" 0 0.1 "10" 0.2 "20", valuelabel) xlabel(,labsize(small)) /// 
			    legend(off)   xtitle("Change in probability (perc. points)") graphregion(color(white)) xtitle(, size(small)) 	///
				order(*.occupation 1.legal_status L   *.residence  *.religion  *.language *.education *.friends  *.team_support *.children 1.country R C P S *.sex) subtitle(A)
				graph save random_effects, replace
				
*combine graphs				
  graph combine  random_effects.gph fixed_effects.gph, graphregion(color(white)) rows(1) xcommon
  graph display 
  graph export "FigS5_fixed-random.png", width(1000) replace	
  graph export "FigS5_fixed-random.eps", replace	
  				
**----------------------------------------------------------------------------------------
**ROBUSTNESS2*
**MM
**no need for interaction terms for margins of non-constrained attributes
   global mod2 = "i.sex  i.religion  i.residence i.occupation i.education i.language i.friends i.team_support i.children  i.country i.legal_status [pweight=weight]"
   qui: reg outcome $mod2, $ses
   margins sex  religion  residence occupation education language friends team_support children, coeflegend post 
   est store m1
*difference in MM country, drop irregular status because of restriction with Romania
	  
       reg outcome i.country if legal_status==2 [pweight=weight], $ses  
	   margins country, coeflegend  post
	   est store m2
*difference in MM legal status, drop romania because of restriction with irregular
       drop if country==3	  
       reg outcome i.legal_status  [pweight=weight], $ses  
	   margins legal_status, coeflegend  post
	   est store m3
*graph
coefplot       (m1, offset(0) msymbol(s) mfcolor(white) color(purple)  ciopts(color(purple)) ) ///	          
			   (m2, offset(0) msymbol(s) mfcolor(white)  color(purple)  ciopts(color(purple))) ///			   			   
			   (m3,  offset(0) msymbol(s)  msymbol(s) mfcolor(white)  color(purple) ciopts(color(purple))), ///
			    headings(1.sex = "{bf: Child gender}" 1.residence = "{bf: Parents' length of residence}" 1.occupation = "{bf: Parents' occupation}"  ///
1.education = "{bf: Parents' education}" 1.religion = "{bf: Parents' religion}" 1.language = "{bf: Parents' Italian proficiency}" 1.friends = "{bf: Family friends Origin}" 1.children = "{bf:Total number of children}" 1.team_support = "{bf:Parents' team support}" 1.country = "{bf: Parents' country of origin}" 1.legal_status = "{bf:Parents' legal status}") ///
			   xtitle(Marginal Means) legend(off) ysize(10) xsize(8) xtitle("Marginal Mean (percentage)", size(small)) ///
			   xlabel(0 0.1 "10" 0.2 "20" 0.3 "30" 0.4 "40" 0.5 "50" 0.6 "60" 0.7 "70" 0.8 "80", valuelabel)xlabel(,labsize(vsmall)) ///
				order(*.occupation 1.legal_status 2.legal_status  *.residence  *.religion  *.language *.education *.friends  *.team_support *.children 1.country 2.country 3.country 4.country 5.country *.sex) ///
			   coeflabels(,labsize(vsmall)) graphregion(margin(zero)) graphregion(color(white))
			   graph export "FigS4_mmaverage.png", width(1000) replace
			   graph export "FigS4_mmaverage.eps", replace
			   
**check country effect with both legal statuses, so without Romania (cases already dropped)
clear all
use "new_data"
drop if country==3
     qui: reg outcome i.country [pweight=weight], $ses
     margins country , coeflegend  post
	   est store m4
	   coefplot       (m4, offset(0) msymbol(s) mfcolor(white) color(purple)  ciopts(color(purple)) )
**----------------------------------------------------------------------------------------

**AMCEs by choice task(1-5)
**AMCEs for first screen
    clear all
    use "new_data"
    keep if screen=="1"
*run analysis
   do "2-AMCEs.do"
*graph
coefplot (reg, keep( *.sex  *.religion  *.residence *.occupation *.education *.language *.friends  *.team_support *.children 1.country) pstyle(p1)) ///  
(Romania, rename ( (1) = R)  \ ///
country2, rename ( (1) = C)  \ ///
country4, rename ( (1) = P) \ ///
country5, rename ( (1) = S) pstyle(p1)) ///
(reg, keep (1.legal_status ) pstyle(p1)) ///
(legal, rename ( (1) = L) pstyle(p1)), ///
headings(1.sex = "{bf: Child gender}" 1.residence = "{bf: Parents' length of residence}" 1.occupation = "{bf: Parents' occupation}"  ///
1.education = "{bf: Parents' education}" 1.religion = "{bf: Parents' religion}" 1.language = "{bf: Parents' Italian proficiency}" 1.friends = "{bf: Family friends Origin}" 1.children = "{bf:Total number of children}" 1.team_support = "{bf:Parents' team support}" 1.country = "{bf: Parents' country of origin}" 1.legal_status = "{bf:Parents' legal status}") ///
xline (-0.4 -0.3 -0.2 -0.1 0.1 0.2, lstyle(grid)) xline (0, lpattern(dash)) ciopts(recast(rcap)) baselevels ///
coeflabels(    R = "Romania"  ///
               C = "China"  ///
               P = "Pakistan" ///
			   S = "Senegal"  ///
			   L = "Regular", labsize(vsmall)) msymbol(s) mfcolor(white)  offset(0.05) /// 
			  graphregion(margin(zero)) xlabel( -0.2 "-20" -0.1 "-10" 0 0.1 "10" 0.2 "20", valuelabel) xlabel(,labsize(small)) /// 
			    legend(off)   xtitle("Change in probability (perc. points)") graphregion(color(white)) xtitle(, size(small))  ///
				order(*.occupation 1.legal_status L   *.residence  *.religion  *.language *.education *.friends  *.team_support *.children 1.country R C P S *.sex)
				graph save "screen1.gph", replace

**AMCEs for second screen
   clear all
   use "new_data"
   keep if screen=="2"

   *run analysis
   do "2-AMCEs.do"
*graph
coefplot (reg, keep( *.sex  *.religion  *.residence *.occupation *.education *.language *.friends  *.team_support *.children 1.country) pstyle(p1)) ///  
(Romania, rename ( (1) = R)  \ ///
country2, rename ( (1) = C)  \ ///
country4, rename ( (1) = P) \ ///
country5, rename ( (1) = S) pstyle(p1)) ///
(reg, keep (1.legal_status ) pstyle(p1)) ///
(legal, rename ( (1) = L) pstyle(p1)), ///
headings(1.sex = "{bf: Child gender}" 1.residence = "{bf: Parents' length of residence}" 1.occupation = "{bf: Parents' occupation}"  ///
1.education = "{bf: Parents' education}" 1.religion = "{bf: Parents' religion}" 1.language = "{bf: Parents' Italian proficiency}" 1.friends = "{bf: Family friends Origin}" 1.children = "{bf:Total number of children}" 1.team_support = "{bf:Parents' team support}" 1.country = "{bf: Parents' country of origin}" 1.legal_status = "{bf:Parents' legal status}") ///
xline (-0.4 -0.3 -0.2 -0.1 0.1 0.2, lstyle(grid)) xline (0, lpattern(dash)) ciopts(recast(rcap)) baselevels ///
coeflabels(    R = "Romania"  ///
               C = "China"  ///
               P = "Pakistan" ///
			   S = "Senegal"  ///
			   L = "Regular", labsize(vsmall)) msymbol(s) mfcolor(white)  offset(0.05) /// 
			  graphregion(margin(zero)) xlabel( -0.2 "-20" -0.1 "-10" 0 0.1 "10" 0.2 "20", valuelabel) xlabel(,labsize(small)) plotregion(margin(l=35)) /// 
			    legend(off)   xtitle("Change in probability (perc. points)") graphregion(color(white)) xtitle(, size(small))  yscale(off)	 ///
				order(*.occupation 1.legal_status L   *.residence  *.religion  *.language *.education *.friends  *.team_support *.children 1.country R C P S *.sex)
				graph save "screen2.gph", replace								

**AMCEs for third screen
   clear all
   use "new_data"
   keep if screen=="3"
   *run analysis
   do "2-AMCEs.do"
*graph
coefplot (reg, keep( *.sex  *.religion  *.residence *.occupation *.education *.language *.friends  *.team_support *.children 1.country) pstyle(p1)) ///  
(Romania, rename ( (1) = R)  \ ///
country2, rename ( (1) = C)  \ ///
country4, rename ( (1) = P) \ ///
country5, rename ( (1) = S) pstyle(p1)) ///
(reg, keep (1.legal_status ) pstyle(p1)) ///
(legal, rename ( (1) = L) pstyle(p1)), ///
headings(1.sex = "{bf: Child gender}" 1.residence = "{bf: Parents' length of residence}" 1.occupation = "{bf: Parents' occupation}"  ///
1.education = "{bf: Parents' education}" 1.religion = "{bf: Parents' religion}" 1.language = "{bf: Parents' Italian proficiency}" 1.friends = "{bf: Family friends Origin}" 1.children = "{bf:Total number of children}" 1.team_support = "{bf:Parents' team support}" 1.country = "{bf: Parents' country of origin}" 1.legal_status = "{bf:Parents' legal status}") ///
xline (-0.4 -0.3 -0.2 -0.1 0.1 0.2, lstyle(grid)) xline (0, lpattern(dash)) ciopts(recast(rcap)) baselevels ///
coeflabels(    R = "Romania"  ///
               C = "China"  ///
               P = "Pakistan" ///
			   S = "Senegal"  ///
			   L = "Regular", labsize(vsmall)) msymbol(s) mfcolor(white)  offset(0.05) /// 
			  graphregion(margin(zero)) xlabel( -0.2 "-20" -0.1 "-10" 0 0.1 "10" 0.2 "20", valuelabel) xlabel(,labsize(small)) /// 
			    legend(off)   xtitle("Change in probability (perc. points)") graphregion(color(white)) xtitle(, size(small)) 	yscale(off) plotregion(margin(l=35))  ///
				order(*.occupation 1.legal_status L   *.residence  *.religion  *.language *.education *.friends  *.team_support *.children 1.country R C P S *.sex)
				graph save "screen3.gph", replace

**AMCEs for fourth screen
   clear all
   use "new_data"
   keep if screen=="4" 

   *run analysis
   do "2-AMCEs.do"
*graph
coefplot (reg, keep( *.sex  *.religion  *.residence *.occupation *.education *.language *.friends  *.team_support *.children 1.country) pstyle(p1)) ///  
(Romania, rename ( (1) = R)  \ ///
country2, rename ( (1) = C)  \ ///
country4, rename ( (1) = P) \ ///
country5, rename ( (1) = S) pstyle(p1)) ///
(reg, keep (1.legal_status ) pstyle(p1)) ///
(legal, rename ( (1) = L) pstyle(p1)), ///
headings(1.sex = "{bf: Child gender}" 1.residence = "{bf: Parents' length of residence}" 1.occupation = "{bf: Parents' occupation}"  ///
1.education = "{bf: Parents' education}" 1.religion = "{bf: Parents' religion}" 1.language = "{bf: Parents' Italian proficiency}" 1.friends = "{bf: Family friends Origin}" 1.children = "{bf:Total number of children}" 1.team_support = "{bf:Parents' team support}" 1.country = "{bf: Parents' country of origin}" 1.legal_status = "{bf:Parents' legal status}") ///
xline (-0.4 -0.3 -0.2 -0.1 0.1 0.2, lstyle(grid)) xline (0, lpattern(dash)) ciopts(recast(rcap)) baselevels ///
coeflabels(    R = "Romania"  ///
               C = "China"  ///
               P = "Pakistan" ///
			   S = "Senegal"  ///
			   L = "Regular", labsize(vsmall)) msymbol(s) mfcolor(white)  offset(0.05) /// 
			  graphregion(margin(zero)) xlabel( -0.2 "-20" -0.1 "-10" 0 0.1 "10" 0.2 "20", valuelabel) xlabel(,labsize(small)) /// 
			    legend(off)   xtitle("Change in probability (perc. points)") graphregion(color(white)) xtitle(, size(small)) ///
				order(*.occupation 1.legal_status L   *.residence  *.religion  *.language *.education *.friends  *.team_support *.children 1.country R C P S *.sex)
				graph save "screen4.gph", replace

**AMCEs for fifth screen
   clear all
   use "new_data"
   keep if screen=="5"

   *run analysis
   do "2-AMCEs.do"
*graph
coefplot (reg, keep( *.sex  *.religion  *.residence *.occupation *.education *.language *.friends  *.team_support *.children 1.country) pstyle(p1)) ///  
(Romania, rename ( (1) = R)  \ ///
country2, rename ( (1) = C)  \ ///
country4, rename ( (1) = P) \ ///
country5, rename ( (1) = S) pstyle(p1)) ///
(reg, keep (1.legal_status ) pstyle(p1)) ///
(legal, rename ( (1) = L) pstyle(p1)), ///
headings(1.sex = "{bf: Child gender}" 1.residence = "{bf: Parents' length of residence}" 1.occupation = "{bf: Parents' occupation}"  ///
1.education = "{bf: Parents' education}" 1.religion = "{bf: Parents' religion}" 1.language = "{bf: Parents' Italian proficiency}" 1.friends = "{bf: Family friends Origin}" 1.children = "{bf:Total number of children}" 1.team_support = "{bf:Parents' team support}" 1.country = "{bf: Parents' country of origin}" 1.legal_status = "{bf:Parents' legal status}") ///
xline (-0.4 -0.3 -0.2 -0.1 0.1 0.2, lstyle(grid)) xline (0, lpattern(dash)) ciopts(recast(rcap)) baselevels ///
coeflabels(    R = "Romania"  ///
               C = "China"  ///
               P = "Pakistan" ///
			   S = "Senegal"  ///
			   L = "Regular", labsize(vsmall)) msymbol(s) mfcolor(white)  offset(0.05) /// 
			  graphregion(margin(zero)) xlabel( -0.2 "-20" -0.1 "-10" 0 0.1 "10" 0.2 "20", valuelabel) xlabel(,labsize(small)) /// 
			    legend(off)   xtitle("Change in probability (perc. points)") graphregion(color(white)) xtitle(, size(small)) 	yscale(off)   ///
				order(*.occupation 1.legal_status L   *.residence  *.religion  *.language *.education *.friends  *.team_support *.children 1.country R C P S *.sex)      
				graph save "screen5.gph", replace

*combine graphs		
   graph combine "screen1" "screen2" "screen3", ///
name("firstset", replace) subtitle("First, second and third choice task", ring(0.1)) graphregion(color(white)) rows(1)   
graph combine "screen4" "screen5", ///
name("secondset", replace) rows(1) xcommon graphregion(margin(l=35 r=30)) graphregion(color(white)) subtitle("Fourth and fifth choice-task", ring(0.1))
graph combine firstset secondset , xcommon ///
saving("screens.gph", replace) rows(2)  graphregion(color(white))
graph export "FigS6.png", replace
graph export "FigS6.eps", replace
		

***--------------------------------------------------------------------------------------------------------------------------------------------






















